NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


THESIS 


VORTEX  RIPPLE  MORPHOLOGY  USING  DUNE2D  MODEL 

by 

Stephen  D.  Martin 

March  2001 

Thesis  Advisor: 

Edward  B.  Thornton 

Thesis  Co-Advisor: 

Timothy  P.  Stanton 

Approved  for  public  release;  distribution  is  unlimited 


20010531  056 


REPORT  DOCUMENTATION  PAGE 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including 
the  time  for  reviewing  instruction,  searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and 
completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any 
other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington 
headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite 
1204,  Arlington,  VA  222024302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project 
(0704-0188)  Washington  DC  20503.  _ 

T  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

Mar  2001  Master's  Thesis 

4.  TITLE  AND  SUBTITLE:  Vortex  Ripple  Morphology  using  DUNE2D  Model  5.  FUNDING  NUMBERS 

6.  AUTHOR(S)  Martin,  Stephen  D.  _ 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESSES)  3.  PERFORMING  ORGANIZATION 

Naval  Postgraduate  School  REPORT  NUMBER 

Monterey,  CA  93943-5000  ■ 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESSES)  10.  SPONSORING  /  MONITORING 
Office  of  Naval  Research  AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official 
policy  or  position  of  the  Department  of  Defense  or  the  U.S.  Government. 

12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT  I 12b.  DISTRIBUTION  CODE 

Approved  for  public  release;  distribution  is  unlimited  (mix  case  letters) 

13.  ABSTRACT  (maximum  200  words) 

The  DUNE2D  boundary  layer  model  for  small-scale  morphology  (Andersen,  1999)  is  compared  with  bedform 
morphology  measured  on  the  inner  shelf  in  11m  water  depth  during  the  SHOWEX  experiment  at  Duck,  N.C.  The 
model  consists  of  inter-linked  modules  for  flow,  sediment  transport  and  morphology.  The  flow  module  is  based 
on  solving  the  Reynold’s  averaged  Navier- Stokes  equations  in  the  vertical  plane  with  k-omega  turbulence  closure. 
The  model  has  been  extended  to  accept  a  general  (but  periodic)  bottom  boundary  to  be  able  to  compare  with  field 
data.  Boundary*  layer  velocity  profiles  were  measured  using  a  Bistatic  Coherent  Doppler  Velocity  profiler 
(BCDV).  A  two-axis  scanning  sonar  altimeter  measured  small-scale  morphology  over  a  1  by  1.5  m  area  with  4  cm 
horizontal  and  0.25  cm  vertical  resolutions.  Bottom  maps  of  small-scale  morphology  were  obtained  continuously 
every  20  minutes.  A  relatively  simple  data  sequence  was  selected  for  model  comparison,  during  which  time  the 
wave  forcing  evolved  abruptly  from  Hsi^=0.3m  to  Hsl?=3.0m  (bed  velocity  <  1  m/s),  and  the  bed  evolved  from  no 
motion  (relic)  to  actively  migrating  vortex  ripples.  SHOWEX  bedform  changes  under  low  wave  plus  collinear 
current  conditions  resulted  in  minor  changes  of  the  vortex  ripple  fields.  Bedform  migration  rates  of  the  model 
were  similar  to  the  field  migration  rates.  Like  the  field  data,  the  modeled  data  under  strong  forcing  removed 
smaller  scale  vortex  ripples  and  redistributed  the  sediment  into  a  larger  scale  ripple  with  a  large  portion  of 
sediments  in  suspension  above  the  bed.  Limitations  of  the  model  owing  to  the  2-D  assumption,  periodic 
boundaries  and  monochromatic  wave  forcing  constraints  are  discussed. 


14.  SUBJECT  TERMS  Waves,  Bottom  boundary  layer,  Nearshore  dynamics,  Sediment  transport,  15.  NUMBER  OF 
Morphology.  PAGES  g^ 

16.  PRICE  CODE 

17.  SECURITY  18.  SECURITY  19.  SECURITY  20.  LIMITATION 

CLASSIFICATION  OF  CLASSIFICATION  OF  THIS  CLASSIFICATION  OF  OF  ABSTRACT 

REPORT  PAGE  ABSTRACT 

Unclassified  Unclassified  Unclassified  UL 

NSN  7540-01-280-5500  Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  239-18 


1 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


11 


Approved  for  public  release;  distribution  is  unlimited 


VORTEX  RIPPLE  MORPHOLOGY  USING  DUNE2D  MODEL 

Stephen  D.  Martin 

Lieutenant  Commander,  United  States  Navy 
B.S.,  United  States  Naval  Academy,  1991 
M.S.,  The  Pennsylvania  State  University,  1993 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  METEOROLOGY  AND  OCEANOGRAPHY 


from  the 


NAVAL  POSTGRADUATE  SCHOOL 
March  2001 


Author: 


Approved  by: 

Edward  B.  Thornton,  Thesis  Advisor 


Timothy  P.  Stanton,  Thesis  Co-Advisor 


Roland  W.  Garwood  Jr.,  Chaii^i 
Oceanography  Department 


iii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


IV 


ABSTRACT 


The  DUNE2D  boundary  layer  model  for  small-scale  morphology 
(Andersen,  1999)  is  compared  with  bedform  morphology  measured  on  the  inner  shelf  in 
1  lm  water  depth  during  the  SHOWEX  experiment  at  Duck,  N.C.  The  model  consists  of 
inter-linked  modules  for  flow,  sediment  transport  and  morphology.  The  flow  module  is 
based  on  solving  the  Reynold’s  averaged  Navier-Stokes  equations  in  the  vertical  plane 
with  k-omega  turbulence  closure.  The  model  has  been  extended  to  accept  a  general  (but 
periodic)  bottom  boundary  to  be  able  to  compare  with  field  data.  Boundary  layer 
velocity  profiles  were  measured  using  a  Bistatic  Coherent  Doppler  Velocity  profiler 
(BCDV).  A  two-axis  scanning  sonar  altimeter  measured  small-scale  morphology  over  a 
1  by  1.5  m  area  with  4  cm  horizontal  and  0.25  cm  vertical  resolutions.  Bottom  maps  of 
small-scale  morphology  were  obtained  continuously  every  20  minutes.  A  relatively 
simple  data  sequence  was  selected  for  model  comparison,  during  which  time  the  wave 
forcing  evolved  abruptly  from  Hsjg=0.3m  to  Hsjg=3.0m  (bed  velocity  <  1  m/s),  and  the 
bed  evolved  from  no  motion  (relic)  to  actively  migrating  vortex  ripples.  SHOWEX 
bedform  changes  under  low  wave  plus  collinear  current  conditions  resulted  in  minor 
changes  of  the  vortex  ripple  fields.  Bedform  migration  rates  of  the  model  were  similar  to 
the  field  migration  rates.  Like  the  field  data,  the  modeled  data  under  strong  forcing 
removed  smaller  scale  vortex  ripples  and  redistributed  the  sediment  into  a  larger  scale 
ripple  with  a  large  portion  of  sediments  in  suspension  above  the  bed.  Limitations  of  the 
model  owing  to  the  2-D  assumption,  periodic  boundaries  and  monochromatic  wave 
forcing  constraints  are  discussed. 
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I.  INTRODUCTION 


This  is  a  study  of  ripple  migration,  evolution,  and  erosion  on  the  inner  continental 
shelf  over  time.  Despite  several  decades  of  analytical  and  empirical  research, 
understanding  the  evolution  of  bedform  morphology  remains  a  fairly  young  field.  In  this 
thesis,  the  morphology  problem  is  considered  for  non-breaking  waves  in  the  shallow-to- 
intermediate  water  depth  range  of  the  inner  shelf.  Current  research  emphasis  is  placed  on 
numerical  simulation  of  this  problem.  This  chapter  provides  motivation  for  study, 
research  objectives,  and  background  on  the  bottom  boundary  layer  with  emphasis  on  the 
interactions  of  hydrodynamic  forcing  on  sediment  transport  and  bedforms. 

A.  MOTIVATION  FOR  STUDY 

Naval  operations  have  rapidly  transitioned  from  “blue”,  open  ocean  into  coastal, 
shallow  and  intermediate  waters.  The  operational  oceanographer  is  required  to  implement 
ocean  models  of  the  littoral  region  to  produce  accurate  forecasts  of  ocean  processes. 
Figure  (1.1)  lists  some  Navy  operations,  as  well  as  dynamic  processes  that  are  of  concern. 


Figure  1.1.  Littoral  Naval  Interests,  (after  www.oc.nps.navy.mil/~stanton/miso/) 
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Over  the  past  decade,  naval  research  has  been  highly  focused  on  adapting  to  this 
new  theater  of  operation.  Field  studies,  laboratory  experiments,  and  numerical  modeling 
research  of  these  highly  spatial  and  temporal  varying  processes  are  motivation  to  improve 
physical  understanding  of  the  waves,  currents,  sediment  transport,  and  bottom 
morphology  in  the  littoral  region.  Advances  are  slow  and  expensive.  Hopefully,  the 
knowledge  gained  through  the  combined  observations  and  simulation  of  the  nearshore 
environmental  processes  will  provide  accurate  oceanographic  forecasts  for  required 
operational  applications. 

Common  to  all  dynamic  processes  in  the  littoral  region  are  the  hydrodynamic 
forces  produced  through  wave,  current,  and  tidal  energy.  Motivation  for  this  thesis  is 
based  on  understanding  wave  and  current  interactions  with  the  bottom  boundary  layer 
that  force  sediment  transport  in  the  form  of  bedload  and  suspended  load,  and  result  in 
bedform  evolution  and  migration. 

Wave-formed  ripples  are  the  dominant  bedforms  in  coastal  regions  (Traykovski, 
et  al,  1999).  Over  time,  mobile  beds  can  evolve  in  such  a  way  that  the  overall  bed 
roughness  is  changed.  Furthermore,  migration  of  sediment  affects  form  drag  on  the 
existing  hydrodynamic  processes,  net  sediment  transport,  and  the  bottom  boundary  layer 
turbulence  generation.  These  changes  can  have  significant  affect  on  wave  energy 
dissipation  (Ardhuin,  2000). 

B.  OBJECTIVES 

The  objective  is  to  simulate  the  migration  and  evolution  of  small-scale  bedforms 

measured  in  the  field  using  DUNE2D,  a  numerical  model  of  the  bottom  boundary  layer 

for  small-scale  morphology.  The  model  is  initiated  using  bedform  and  velocity 
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measurements  acquired  during  the  SHOWEX  field  experiment  conducted  from 
September  to  December  1999  at  Duck,  North  Carolina.  Small-scale  bottom  maps  and 
high-resolution  velocity  measurements  of  a  specific  storm  event  are  used  to  qualitatively 
verify  DUNE2D  output.  Efforts  are  focused  on  simulating  actual  field  conditions  from 
an  initial  bed  configuration.  The  long-term  objective  is  to  provide  model  validation  to 
enhance  numerical  prediction  of  nearshore  processes,  particularly  to  predict  ripple 
formation  and  migration. 

C.  BACKGROUND 

1.  Flow  Properties 

Accurately  determining  the  flow  field  in  the  bottom  boundary  layer  is  essential  to 
predicting  sediment  transport  and  morphology  evolution.  Flow  in  the  ocean  is  by  nature 
turbulent.  Sediment  transport  calculations  require  estimates  of  the  bed  stress,  Xb, 

rb=pu;  (1.1) 

where  the  shear  stress  is  proportional  to  the  frictional  velocity,  u*2.  Flow  field  dynamics 
are  complicated  by  irregular  waves.  In  the  inner  shelf  regions,  flow  near  the  bed  is 
influenced  by  wave  groupiness.  Quantitative  representation  is  made  difficult  through 
multiple  wave  packets  interacting  with  each  other  superimposed  with  mean  currents  of 
varying  strength  and  direction.  Deigaard  et  al  (1999)  found  net  sediment  transport  under 
regular  waves  and  wave  groups  with  bound  long  waves  to  be  dependent  on  sediment  size, 
as  well  as  wave  characteristics. 
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2.  Sediment  Transport  through  the  Shields  Parameter 

Sediment  transport  can  be  described  in  terms  of  the  Shields  parameter.  The 
Shields  parameter  is  the  ratio  of  the  destabilizing  forces  of  the  bottom  stress  and  the 
stabilizing  force  of  gravity  acting  on  the  sediment, 
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where  g  is  the  gravitational  acceleration,  d  is  the  sediment  diameter,  and  s  is  the  ratio  of 
sediment  density  to  fluid  density.  Sediments  have  a  critical  shear  stress  that  is  dependent 
upon  their  grain  size,  density,  porosity,  and  shape.  Once  the  shear  stress  exceeds  this 
critical  threshold,  the  viscous  sub-layer  breaks  down,  and  the  sand  grains  are  set  in 
motion.  As  the  strength  of  the  flow  increases,  the  grains  are  lifted  off  the  bed.  As  long  as 
the  settling  velocity  of  the  sediment  is  small  compared  with  the  forces  producing  lift,  the 
grains  will  remain  in  suspension  above  the  bed.  The  grains  redistribute  back  to  the  bed 
when  the  fluid  forces  acting  on  the  grains  become  less  than  the  gravitational  forces. 
Deposition  is  difficult  to  predict,  as  the  flow  field  has  eddies  of  varying  strength  and  size. 
For  high  Shields  numbers  (0>1)  under  non-breaking  waves,  Diegard  et  al  (1999)  found 
that  sheet  flow  occurs,  where  a  large  portion  of  the  bedform  sediment  is  set  into  motion, 
such  that  the  remaining  bed  is  planar  in  shape.  For  forcing  with  low  Shields  numbers, 
ripple-dominated  bedforms  occur. 

3.  Bedforms 

Ripples  are  classified  by  their  geometry  and  commonly  divided  into  two  separate 
spatial  scales.  Dunes,  antidunes,  megaripples  and  relic  ripples  are  considered  large-scale 
bedforms.  Small-scale  ripples  include  vortex,  cross-ripples  and  rolling  grain  ripples. 
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Small-scale  ripples  are  categorized  by  a  length  of  less  than  0.6  meters  and  a  height 
roughly  ten  percent  of  the  length,  or  below  0.06  meters  (Engelund  and  Fredsoe,  1982). 

Bagnold  (1946)  described  two  types  of  wave-induced  ripples:  rolling  grain  and 
vortex  ripples.  For  increasing  Shields  numbers,  as  the  critical  value  is  just  exceeded  for 
flow  over  a  flat,  sandy  seabed,  grains  oscillate  back  and  forth  in  a  rolling  motion.  This 
condition  is  difficult  to  find  in  nature  (Andersen,  1999).  As  the  flow  intensity  increases 
and  the  turbulent  eddies  induce  more  sediment  transport,  the  oscillatory  wave  action 
organizes  the  moving  grain  into  the  classic  vortex  ripple  triangular  shapes.  Over  time, 
the  sediment  is  moved  from  the  trough  of  the  ripple,  and  deposited  through  bedload  and 
suspended  load  towards  the  ripple  peak.  The  steepness  of  the  ripple  increases. 
Eventually,  the  flow  vortex,  formed  by  the  fluid  moving  over  the  ripple  peak,  separates 
from  the  bed  and  induces  further  sediment  transport  (Neilsen,  1981).  Many  studies  have 
been  conducted  to  find  equilibrium  or  quasi-equilibrium  forms  of  the  bed  due  to  wave 
and  current  forcing  with  the  objective  of  predicting  ripple  wavelength,  height,  and 
steepness  as  a  function  of  mobility  number  or  Shields  parameter. 

4.  Numerical  Solutions 

Several  numerical  approaches  have  been  used  to  describe  the  flow  in  the  bottom 
boundary  layer  over  an  arbitrary  bedform.  Perrier’s  (1996)  discrete  vortex  (DV)  method 
does  a  good  job  describing  the  qualitative  features  of  the  flow.  Turbulence  models 
provide  better  quantitative  description,  which  is  required  for  morphological  calculations 
(Andersen,  1999).  The  k-s  turbulence  closure  technique  was  used  by  Sato  (1986), 
Tsujimoto  et  al.  (1991),  and  Perrier  (1996),  while  Wilcox  (1988,  1993b)  introduced  the 
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k-o  turbulence  closure.  The  k-cd  model  handles  the  strong  pressure  gradients  near  the 
separation  point  of  the  turbulent  eddies  better  than  the  k-s  model. 

The  present  DUNE2D  model  study  uses  the  k-o  closure.  The  model  is  presented 
in  the  next  chapter.  Good  results  have  been  obtained  by  Andersen  (1999)  using 
DUNE2D  for  monochromatic  wave  and  wave  plus  current  forcing  compared  with 
laboratory  data  (Fredsoe  et  al,  1999). 
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II.  DUNE2D  MODEL  FORMULATION 


DUNE2D  is  a  bottom  boundary  layer  model  for  sediment  transport  over  small- 
scale  ripple  morphology  (Andersen,  1999).  The  model  resolves  two-dimensional  flow 
over  a  one-dimensional  bedform.  DUNE2D  is  a  modular  Fortran  code,  which  evolved 
over  many  years  of  research  on  flow,  sediment  transport,  and  most  recently,  morphology 
at  the  Technical  University  of  Denmark.  DUNE2D  is  composed  of  three  modules:  flow, 
transport,  and  morphology.  A  flow  diagram  of  the  model  is  represented  in  Figure  (2.1). 


Figure  2. 1 .  DUNE2D  Model  Schematic  Flow  Diagram. 


The  program  has  several  modes  of  simulation.  For  each  mode,  several  physical 

assumptions  and  numerical  modeling  techniques  exist.  The  model  can  be  run  for  flow 

and  sediment  transport  simulations,  or  for  bottom  boundary  layer  flow  alone.  For 

purposes  of  this  study,  systematic  solutions  for  all  three  principle  modules  are  considered 

with  emphasis  on  the  evolution  of  bedforms.  Therefore,  the  model  description  in  this 

chapter  will  be  presented  from  a  morphology  simulation  perspective. 
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Dynamic  forcing  is  through  periodic,  sinusoidal  wave  input.  If  desired,  a  mean 
along-axis  current  can  be  included.  At  each  grid  point  of  the  model,  the  horizontal,  u, 
and  vertical,  w,  velocity  flow  components  are  evaluated.  Flow  solutions  are  passed  to  the 
transport  module  where  sediment  characteristics,  as  well  as  the  flow  solutions,  influence 
sediment  transport  calculations.  Strong  turbulent  flow  is  capable  of  sediment  transport, 
both  along  the  bed  and  in  suspension  above  the  bed.  Transport  solutions  are  passed  to  the 
morphology  module  where  bedform  changes  are  integrated  over  time. 


DUNE2D  models  oscillatory  flow  with  an  optional  mean  collinear  current  over  a 
rippled  bedform,  as  illustrated  by  Figure  (2.2).  The  twelve  parameters  listed  in  Figure 
(2.2)  are  required  to  describe  the  model  domain.  Appendix  (A)  provides  a  further 
explanation  of  the  input  parameters.  Hydrodynamic  forcing,  geophysical  properties  of 
the  sediment,  numerical  schemes,  grid  definitions,  boundary  conditions,  transport 
methods,  and  morphology  setup  are  some  of  the  categories  that  are  used  to  describe 
model  input.  Appendix  (B)  provides  examples  of  the  ASCII  input  files  and  the 
procedural  method  used  to  implement  DUNE2D  code. 


v:  kinematic  viscosity 
p:  fluid  density 
ps:  sediment  density 
X:  ripple  length 
h:  ripple  height 
d^:  sediment  grain  diameter 
D:  depth  of  flow 
U0:  amplitude  of  oscillatory 
wave  velocity7 

T:  period  of  oscillatory  wave 
Uc:  mean  current 
g:  gravitational  acceleration 
k>:  Nikuradse  bed  roughness 


Figure  2.2. 


Flow  over  a  ripple  bed. 
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A.  FLOW  MODULE 

The  flow  module  models  fully  turbulent  flows,  with  Reynolds  numbers  ranging 
between  105  to  106 ,  in  the  bottom  boundary  layer.  The  model  initializes  input  variables 
and  sets  necessary  default  values,  and  then  solves  the  flow  field  given  specific  dynamic 
forcing  over  an  arbitrary  one-dimensional  bedform. 

1.  Conservation  of  Momentum  and  Mass 

The  two-dimensional  flow  equations  along  the  horizontal  and  vertical  axis  assume 
an  incompressible,  Newtonian  fluid.  The  conservation  of  momentum  is, 

dw  dui  l  dp  d  f  1  f  dui  dUj  ^  /2i\ 

dt  3  dXj  p  dxi  dxi  ^2|^<3x;.  dxi  J 

where  p,  uh  p,  and  v  represent  pressure,  velocity,  density  and  kinematic  viscosity, 
respectively.  The  continuity  equation  is  described  by, 

#  =  0  i  =  1,2  (2.2) 

ox, 

The  Reynolds-averaged  Navier-Stokes  (RNS)  momentum  equations  are  obtained 
by  applying  Reynolds  decomposition  and  averaging  to  Equations  (2.1)  and  (2.2)  with  the 
pressure  and  flow  quantities  divided  into  a  mean  and  fluctuating  components  (e.g. 
I4j  Ui+Uj  ). 
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where  the  mean  component  resolves  the  waves  and  currents,  denoted  by  a  capital  letter, 
the  turbulence  fluctuations  are  denoted  by  primes,  and  the  over  bar  indicates  time 
averaging.  A  slow  time  variation  of  the  mean  velocity  is  allowed. 

The  last  term  in  Equation  (2.3)  is  the  Reynolds  stress  tensor.  The  RNS  equations 
are  closed  using  turbulent  viscosity  closure,  where  the  stress  tensor  is  modeled  as 


f'fdU.  dU,  ^ 
!  +  ; 


-UjUi=VT 


K-KdXj  dx, 


JJ 


dx 


-~kSif 
3  * 


(2.5) 


where  k  =  —  (w  2  +w'2 ) 

vT  is  the  turbulent  viscosity  and  6y  is  the  Kronecker  delta.  Equation  (2.5)  describes  a 
dampening  and  diffusion  process. 

2.  Turbulent  Kinetic  Energy  Equations  using  k-©  Model 


A  two-equation  k-q  model  is  used  to  determine  eddy  viscosity  as  a  function  of 
time  and  space  (Anderson,  1999).  The  k-cd  model  (Wilcox  (1988;  1993b))  was  selected 
over  the  more  commonly  used  k-s  model  to  accommodate  boundary  conditions  on  a 
rough  bed  in  unsteady  flow  with  strong  pressure  gradients  along  the  bottom  boundary. 
The  turbulent  viscosity  term  in  Equation  (2.5)  has  a  direct  relationship  to  k  (kinetic 
energy)  and  a  (vorticity)  through  a  model  constant  y  . 


* 


vT=r 


K 

VJ 


(2.6) 
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The  turbulent  kinetic  energy  equations  are  used  to  solve  for  k  and  a.  (Wilcox,  1993b) 


8k  5k  _  d 
dt  1  dx]  dXj 


(v  +  ovT) 


Sk 

dx. 


- dUi  /?* 

~U‘Uj KXS 


(2.7) 


dm  TT  dm  d 

—  +  U,- —  =  — 

dt  ;  6xj  dXj 
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-  firn1 


(2.8) 


where  closure  coefficients :  y*  =  1;  y  -  — ;  0*  =  — — ; 

9  100 

n  3  1,1 

H  40  2  2 

These  equations  are  similar  to  the  k-s  model,  where  s  is  the  local  dissipation  rate 
of  turbulent  kinetic  energy.  The  term  ©  does  not  have  a  direct  physical  meaning,  but  can 
be  described  as  the  vorticity  of  turbulent  eddies  (Saffmen,  1970),  and  is  related  to  e  by, 


s  =  0*  Km 


(2.9) 


3.  Roughness  of  the  Wall 

The  bottom  boundary  is  scaled  using  the  concept  of  a  rough  wall  through  the 
Nikuradse  roughness,  kN/30.  In  DUNE2D,  kn  is  set  to  two  and  a  half  times  the  median 
grain  diameter.  (Fredsoe  and  Diegaard,  1992). 

kN  =  2.5d50  (210> 

The  law  of  the  wall  is  applied  to  the  bottom  boundary  condition  associated  with 
the  sandy  bed.  A  modification  is  made  to  a  smooth  wall  to  fit  the  empirically  rough  wall 
solution.  (Andersen,  1999).  The  bottom  boundary  condition  is  a  no-slip  condition  for  the 
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flow,  which  implies  a  viscous  sub-layer  and  k=0  at  the  bed,  while  ©  follows  (Wilcox, 
1993b): 


(2.11) 


(2.12) 


Andersen  states  that  the  methodology  for  attaining  a  rough  wall  solution  through 
parameterization  constant,  Sr,  is  not  physically  correct.  Wall  functions  are  normally  used 
with  turbulent  kinetic  energy  models  that  don’t  involve  viscous  sub-layers.  The  k-co 
model  has  a  viscous  sub-layer  due  the  no-slip  condition  imposed  at  the  bottom  boundary. 
Andersen  also  points  out  that  model  results  using  the  k-co  model  are  sensitive  to  the  free- 
stream  boundary  condition  for  co.  To  overcome  this  problem,  a  symmetric  boundary  is 
applied  to  the  upper  boundary.  (Andersen,  1999) 

B.  SEDIMENT  TRANSPORT 

DUNE2D  separates  sediment  transport  into  bedload  and  suspended  load 

calculations  using  techniques  described  by  Engelund  and  Fredsoe  (1976)  and  Fredsoe  and 

Diegaard  (1992).  Turbulent  flow  causes  destabilizing  forces  upon  the  rippled  bed. 

Sediment  is  thrown  into  motion  when  strong  destabilizing  forces,  specifically,  drag  and 

lift,  overcome  gravitational  stabilizing  forces.  Segregation  of  bedload  and  suspended 

load  is  imposed  by  an  imaginary  boundary.  Moving  sediment  in  repeated  contact  with 
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the  bed  is  known  as  bedload  transport.  Suspended  load  transport  is  sediment  in  motion 
above  the  bedload  transport  boundary,  defined  at  a  height  2d,  by  a  function 


z=2d 


(2.13) 


where  0  is  the  Shields  parameter  and  y  is  the  slope  of  the  bed  and  d  is  the  median  grain 
diameter.  Both  transport  load  calculations  are  based  on  the  Shields  parameter,  0 ,  for  a 
flat  bed  (Equation  1.2),  where  the  maximum  shear  stress  during  a  wave  cycle  is 


where  U0  is  the  wave  velocity  amplitude  and  the  wave  frictional  coefficient  is, 


(2.14) 


f  a  V0  25 


L  =  0.04 


\Kj 
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where  A  is  the  orbital  amplitude  of  the  wave  motion. 


(2.15) 


A..uj  (2.16) 

2  7t 

The  local  bed  shear  must  be  greater  than  the  critical  bed  shear  for  sediment  motion  to  be 
initiated.  Bedload  is  a  function  of  the  bed  shear,  while  suspended  load  is  a  function  of 
bed  shear  as  well  as  settling  velocity.  The  settling  velocity  of  the  sediment  is  determined 
using  an  empirical  formula  for  the  drag  coefficient  of  sand  (Fredsoe  and  Deigaard,  1992): 

(2.17) 


CD  =1-4  + 
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where  Cd  is  the  drag  coefficient  and  Red  is  the  sediment  Reynolds  number. 

1.  Bedload  Calculations 

Bedload  calculations  are  based  on  the  Engelund  and  Fredsoe  (1976)  formulation. 
Sediment  velocity  is  found  from  the  balance  of  the  gravity,  drag,  and  friction  forces  on  a 
single  grain  on  a  sloping  bed.  The  non-dimensional  single  grain  velocity  is  given  by, 

;  a  - 10  (emperical  constant);  (2.18) 


Ei 

u, 


—  -a 


r  | 

0  \ 

1-0.7  I 
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e 

where  Ub  is  the  grain  velocity  and  9cy  is  the  critical  Shields  number  for  a  sloping  bed, 


Ocy  =  Oc 


f  ■  \ 

sinr 

a  cosy  + - 

V  ) 


(2.19) 


in  which  and  o  is  the  sign  of  the  direction  of  the  bed  load  and  the  critical  Shields 


parameter  for  a  flat  bed  is. 


0  =  4Md 
c  3  CDa2 


for  which  pD,  is  the  dynamic  friction  coefficient.  (Andersen,  1999) 

It  is  assumed  only  a  single  layer  of  grains  can  be  in  motion  at  the  bed.  From  this 
assumption,  an  equation,  based  on  empirical  data,  is  used  to  solve  the  fraction  of  grains 
that  are  in  motion  per  unit  area. 


- 

4 

i+ 

7C 

o  -Qc 

l  J 

(2.20) 
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Knowing  the  number,  or  probability  of  grains  in  motion,  bedload  transport,  or  bedload 
flux,  the  rate  of  transport  of  a  volume  of  sediment  per  time,  is  determined.  The  non- 
dimensional  bedload  flux  is  (Fredsoe  and  Diegaard,  1992) 


After  substituting  for  n,  Ub  ,  and  including  Fredsoe’ s  linear  gravity  correction  to  the 
Meyer-Peter  flat  bed  equation  for  bedload  flux,  Equation  (2.21)  simplifies  to: 

h(0,r)=s(e-0'-Prf!  (2'22) 

where  (3  =  0.1. 

2.  Suspended  Load  Calculations 

Suspended  sediment  is  advected  by  wave  motion  and  mean  currents  while  falling 
towards  the  bed  based  on  the  sediment’s  settling  velocity,  ws.  Strong  turbulent 
fluctuations  with  small  settling  velocity  sediment  produce  sediment  in  suspension  above 
the  bedload  boundary.  Suspended  sediment  concentration,  c,  is  modeled  by, 

iY3 

c  =  c0  1+—  ;  c0  =  0.65 gm  >  (2.23) 

l  Af) 

where  \f,  the  mean  free  path  of  the  sediment  at  the  bedload  boundary  condition,  is 


4  y2 

0.013  ps6 


(2.24) 


and  n,  again,  represents  the  fraction  of  grains  in  motion  per  unit  area  as  described  by 
Equation  (2.20). 
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The  total  instantaneous  suspended  concentration  is  found  using  Equation  (2.23) 
where  the  time  rate  change  of  concentration  is  calculated  using  a  transport-diffusion 

equation: 

Dc  dc  _  n  (2.25) 

—  =  w—+V{stVc) 

Dt  oz 

where  the  settling  velocity  is  represented  by  wt  using  Equation  (2.17)  and  the  term  es  *vT 
is  the  suspended  sediment  diffusivity,  which  is  equal  to  the  eddy  viscosity  defined  by 
Equation  (2.6). 

The  instantaneous  suspended  sediment  flux  is  obtained  through  the  instantaneous 
velocity  distribution  and  concentration  found  from  Equation  (2.25), 

qs  =  u(x,  z,  f)c(x,  z,  t )  •  (2-26) 


Finally,  The  total  sediment  load  flux,  qt,  through  a  cross  section,  is  the  vertical  integral  of 
the  flux  profile. 


q,(x,t)  qb(x,t)+ 


J q,(x,z,t)dz 

h(x,t) 


(2.27) 


Andersen  states  that  one  of  the  fundamental  problems  with  the  calculation  of  the  bedload 
transport  is  that  it  is  based  on  mean  bed  shear  stresses.  The  model  does  not  incorporate 
the  turbulent  fluctuations,  k.  Despite  this  deficiency,  dynamic  ripple  mechanisms  are 
preserved.  (Andersen,  1999) 
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C.  MORPHOLOGY  MODULE 


The  morphology  module  is  based  on  the  continuity  of  mass  for  sediment.  From 
the  total  transport  Equation  (2.27),  the  change  of  the  bedform  over  time  is  found  by  time 
integrating, 

(2.28) 

dhjxj)  _ _ 1  dgXx,  t) 

dt  1  -  p  dx 

where  p  is  the  porosity  of  the  sand.  The  main  interest  for  the  morphology  simulations  is 
the  mean  bed  height  at  the  end  of  each  wave  period.  For  further  discussion  of  the  output 
parameters  see  Appendix  (A). 

D.  WAVE  PLUS  CURRENTS  SUPPLEMENT 

The  model  forcing  can  be  supplemented  with  the  addition  of  a  collinear  current  on 
the  oscillatory  flow.  The  addition  of  the  current  is  made  through  a  replacement  of  the 
constant  friction  factor  for  wave-only  flows  to  a  wave  plus  current  friction  factor,^;;, 


cup  Y 

Q  J 


(2.29) 


where /wc  is  the  constant  friction  factor  for  the  current  in  a  wave  plus  current  flow.  D  is 
the  depth  scaling  term  and  Q  is  the  water  discharge  per  unit  width  (Fredsoe  et  al.,  1999). 
An  equivalent  Nikuradse  roughness,  kwc,  is  related  to  fwc  by 


f 

14.8Z)exp 

V 


■0.57 


(2.30) 


The  advection  of  sediment  by  the  mean  flow  is  calculated  in  the  bedload  and  suspended 
load  modules. 
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III.  FIELD  DATA 


A  long-term  inner  shelf  study  of  the  bottom  boundary  layer  was  part  of  the 
Shoaling  Waves  Experiment,  SHOWEX,  at  Duck,  North  Carolina  (Figure  3.1).  Multi¬ 
scale  measurements  of  the  hydrodynamic  interactions  over  sandy  rippled  bedforms  were 
made.  This  chapter  describes  the  data  collection  procedure,  measurement  techniques,  and 
bedform  morphology  due  to  ocean  forcing  caused  by  a  specific  storm  event. 


Figure  3.1.  U.S.  Army  Corps  of  Engineers  Field  Research  Facility  at  Duck,  North 

Carolina.  Site  of  SHOWEX  99.  (from  http://www.fff.usace.army.mil/) 

A.  MEASUREMENTS 


SHOWEX  was  an  Office  of  Naval  Research  funded  Directed  Research  Initiative 

(DRI)  to  improve  wave  propagation  models  in  the  coastal  regions.  One  of  the  objectives 

was  to  measure  the  dissipation  of  shoaling  waves  across  the  inner  shelf  with  emphasis  on 

the  affects  of  the  feedback  between  wave  induced  ripples  and  the  resulting  increased 
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dissipation  of  the  waves  at  the  bed.  Measurements  of  the  bottom  boundary  layer  were 
made  at  11  meters  water  depth  from  a  bottom  mounted  platform  with  precision 
instruments  at  site  S-3  (Figure  3.2)  at  Duck,  North  Carolina,  during  28  September  to  10 
December  1999. 


y$,* 


Figure  3.2.  SHOWEX  99  measurement  platform  location  along  Duck,  N.C.  (after 
http://cheyenne.rsmas.miami.edu/duck99/ )  The  platform  is  approximately  1.4km  off  the 
beach. 

Sediment  collection  of  the  relatively  homogeneous  fine  sand  was  made  in  the 
observation  area  four  times  during  the  three-month  experiment.  The  median  sediment 
diameter  was  0. 1  mm  based  on  a  sieving  procedure. 

A  schematic  representation  of  the  instrument  platform  is  shown  in  Figure  (3.3), 
illustrating  the  sensors  used  to  measure  velocities  within  the  bottom  boundary  layer  as 
well  as  the  bedform  and  suspended  sediment  concentrations.  A  1 .3  MHz  pulsed,  bistatic 
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Figure  3.3.  Bottom  boundary  layer  measurement  platform. 


coherent  Doppler  velocity  (BCDV)  profiler  (Stanton  1996,  2000)  measured  vertical 
profiles  of  three  component  velocities,  u ,  v,  and  w,  and  sediment  concentration  in  2  cm 
diameter  volumes  of  water  every  0.6  cm  over  the  bottom  60  cm,  starting  at  30  cm  from 
the  sensor  face.  The  sampling  frequency  was  1 8Hz,  resulting  in  high  resolution  velocity 
profiles.  The  BCDV’s  orientation  is  illustrated  in  Figure  (3.4).  Tilt  sensors  in  the  BCDV 
allowed  precise  orientation  of  the  vertical  coordinate  system.  The  raw  BCDV  vector 
velocity  profile  data  were  transformed  into  an  earth-referenced  coordinate  system  during 


post-processing. 
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F’gure  3.4.  Platform  measurement  orientation  is  22°  to  the  left  of  Beach  North. 
Platform  located  in  1  lm  depth. 

A  two-axis,  scanning,  high-resolution  acoustic  altimeter  was  collocated  with  the 
BCDV  profiler  on  the  instrument  platform.  The  scanning  altimeter  mapped  a  4m  by  2m 
bed  area  every  15  to  30  minutes.  Resolution  of  the  bed  was  4cm  in  the  horizontal  and 
0.25  cm  in  the  vertical  over  an  inner  1  m2  area.  Horizontal  resolution  decreases  away 
from  the  measurement  center  due  to  step  sensor  angles  and  distances  of  travel.  Post¬ 
processing  of  the  raw  range  data  included  removal  of  false  targets,  transformation  of 
range  and  tilt  angles  into  Cartesian  coordinates,  and  objective  analysis  techniques  to 
smoothly  map  the  scanned  data  points  across  the  measurement  area.  Final  bedform 
observation  boundary  limits,  shown  in  Figure  (3.5),  were  reduced  to  1.5m  by  lm  box 
defining  an  area  of  maximum  resolution. 
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SHOWEX  ‘99.  Black  line  is  the  BCDV  measurement  center.  Right  foreground  is  a  hole 
generated  by  a  live  shellfish  embedded  in  the  ocean  floor. 

B.  VELOCITY  FORCING  AND  BEDFORM  MORPHOLOGY  DURING  11-13 

NOVEMBER  1999 

The  evolution  and  migration  of  bedform  are  examined  for  a  storm  event  during 
year  days  315.5  to  317  (November  11-13).  The  storm  produced  strong  local  winds  that 
caused  the  significant  wave  height  to  increase  from  0.3m  around  year  day  315.5  to  a  peak 
of  nearly  3  m  at  year  day  316.4,  as  measured  by  the  FRF  8m  bottom  pressure  gauge  array 
(Figure  3.6).  During  the  storm,  wave  forcing  increased  the  orbital  displacement  just 
above  the  bottom  boundary  layer  from  0.1m  to  1.1m.  As  the  storm  tracked  out  of  the 
region,  local  wind  waves  and  swell  dissipated. 
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data  collection  for  November  1999.  Box  indicates  the  storm  event  of  10-13  November, 
(after  http://www.frf.usace.army.mil/) 


The  u  and  v  velocity  components  measured  by  the  BCDV  were  arithmetically 
averaged  over  the  vertical  between  10cm  to  30cm  above  the  bed  to  form  278  second  time 
series  at  0.1  year  day  intervals  (2.4  hours).  The  velocities  represent  the  near-bed 
velocities  just  above  the  turbulent  boundary  layer,  sometimes  referred  to  as  ux  and  vx. 
From  the  u  and  v  time  series,  the  wave  oscillatory  velocity  amplitude,  U0,  peak  wave 
period  and  mean  currents  were  determined  (Figure  3.7). 
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Figure  3.7.  Mean  currents  and  oscillatory  wave  velocity  amplitude  for  year  days  315.5 
to  3 16.6  illustrating  storm  event  spin-up  of  the  ocean  dynamics. 


The  wave  velocity  amplitude  is  calculated  as  the  equivalent  sinusoidal  velocity 
amplitude  that  gives  the  same  wave  variance  as  measured, 

u.  =  (3-D  . 

where  ou2  and  cv2  are  the  along-shore  and  cross-shore  velocity  component  variances. 
The  peak  wave  period  of  the  initially  mild  waves,  composed  of  local  and  distant  swell 
was  not  well  defined.  By  year  day  315.9  the  storm  peak  wave  period  became  well 
defined  and  increased  over  the  following  24  hours.  The  majority  of  the  mean  current 
energy  was  in  the  along-shore  direction  as  the  magnitude  increased  from  15  cm/s  to  a 
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peak  of  45  cm/s  down-coast.  The  cross-shore  current  remained  relatively  steady  at  10 
cm/s  in  the  onshore  direction. 


Bed  changes  over  time  obtained  from  a  centerline  cross-shore  transect  of  the 
scanned  acoustic  altimeter  data  are  shown  in  Figure  (3.8).  Initially  at  year  day  315.5,  the 
two-dimensional,  cross-shore  transect  bedform  can  be  described  as  small-scale  vortex 
ripples  overtop  a  larger  scale  ripple.  The  bed  became  mobile  as  the  dynamic  forcing 
increased.  By  the  peak  of  the  wave  and  current  forcing,  the  small-scale  ripples  have  been 
planed-off  and  the  larger  scale  feature  has  been  flattened. 
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Figure  3.8.  Evolution  of  bedform  along  centerline,  cross-shore  transect  during  storm 
event.  (Positive  x  direction  is  towards  the  beach). 
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IV.  DUNE2D  MORPHOLOGY  OUTPUT  ANALYSIS 


Sensitivity  tests  of  DUNE2D  morphology  simulations  were  first  performed  to 
determine  model  output  tendencies  given  specific  monochromatic  wave  forcing  and  wave 
plus  current  forcing.  The  model  morphology  evolution  and  migration  is  then  compared 
with  field  data  observed  during  a  storm  event  of  SHOWEX  ‘99.  In  both  studies,  a  cross¬ 
shore  transect  bedform  was  chosen  from  the  SHOWEX  ’99  altimeter  data  and 
transformed  into  a  periodic  boundary  bedform  represented  in  Figure  (4.1).  All 
simulations  use  fine  sediment  grain  diameter  equal  to  0.1mm,  as  measured  at  SHOWEX. 


Figure  4.1.  Transformed  cross-shore  transect  of  bedform  from  SHOWEX  ’99  data 
(year  day  315.9)  used  to  initialize  DUNE2D  simulations.  The  beach  is  to  the  right  with 
the  initial  oscillatory  wave  forcing  going  towards  the  right  (See  Appendix  (C)  for 
explanation  of  bedform  transformation). 


A.  MODEL  SENSITIVITY  TESTS 

Sensitivity  studies  of  DUNE2D  morphology  simulations  are  conducted  on  the 
initial  bedform  (Figure  4.1)  using  varying  initial  wave  forcing  and  wave  plus  current 
collinear  forcing  over  a  short  simulation  time  duration  of  20  wave  periods  to  obtain  a 
qualitative  feeling  for  the  model  performance. 
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1.  Wave-only  Case  Study 


Four  wave-only  forcing  tests  were  made  with  different  monochromatic  wave 
input  values  as  summarized  by  Table  (4.1).  Tests  1  and  2  represent  the  lower  wave 
Table  4. 1 ,  Four  test  cases  in  wave-only  sensitivity  study  of  morphology  simulations. 


Wave  Input 

Test  #1 

Test  #2 

Test  #3 

Test  #4 

Oscillatory  Amplitude,  U0 
(m/s) 

0.35 

0.35 

0.70 

0.70 

Period,  T  (s) 

4.8 

10.0 

4.8 

10.0 

Orbital  Displacement  A,  (m) 

0.23 

0.56 

0.53 

1.1 

Mobility  Number,  \| / 

76 

76 

303 

303 

Field  Measured  Bedform 
(Year  day) 

315.9 

315.9 

315.9 

315.9 

energy  simulations,  with  a  mobility  number,  \|/,  equal  to  76,  indicative  of  the  pre-storm 
condition.  Tests  3  and  4  are  the  higher  wave  energy  simulations  occurring  during  the 
storm  event,  with  a  \|/=303. 

Bedform  changes  over  twenty  wave  periods  for  varying  monochromatic 
oscillatory  wave-only  forcing  is  shown  in  Figures  (4.2a-d).  In  each  figure,  the  output 
bedforms  for  the  first  and  twentieth  period  are  displayed.  Similar  bedform  changes 
resulted  for  test  cases  1  and  2,  while  test  cases  3  and  4  are  similar.  Vortex  ripples  are 
created  over  top  the  larger  scale  bedform  in  test  cases  1  and  2.  The  amplitude  of  the 
vortex  ripples  is  greater  for  the  longer  period  waves  of  case  2.  The  wavelengths  are  the 
same  for  both  cases  with  the  average  distance  between  the  new  ripple  peaks  is  equal  to 
10cm. 
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Figure  4.2.  DUNE2D  bedform  morphology  for  field  bedform  year  day  315.9  with 
wave-only  monochromatic  forcing  after  20  wave  periods.  (a)Test#l:  Uo=0.35m/s, 
T=4.79s,  \(/= 76;  (b)  Test#2:  U0=0.35  m/s,  T=10s,  y=76;  (c)  Test#3:  Uo=0.70m/s, 
T=4.79s,  \|/=303;  (d)  Test#4:  U0=0.70m/s,  T=10s,  \j/=303. 


For  the  high  mobility  number  cases  3  and  4,  the  vortex  ripples  are  planed  off, 
leaving  a  smoother  large-scale  bedform.  The  model  tends  to  quickly  change  the  bed  after 
a  short  time  equal  to  5-10  wave  periods. 

In  summary,  low  energy  wave  forcing  used  in  tests  1  and  2  results  in  the  creation 
of  vortex  ripples  on  top  of  the  larger  scale  bedform.  The  high  energy  wave  forcing  used 
in  tests  3  and  4  planes  off  small-scale  ripples  on  top  of  the  larger  scale  ripple  bedform. 

2.  Wave  plus  Collinear  Current  Case  Study 

The  effects  of  adding  a  collinear  current  to  the  same  low  amplitude,  short  period 
wave  forcing  and  initial  field  measured  bedform  (SHOWEX  ’99  year  day  315.9)  are 
examined.  Four  separate  collinear  currents  (Table  4.2)  are  considered  where  the  ratio  of 
current  to  wave  velocity  amplitude  is  Uc/U0  —  0.5  (Test  5),  Uc/U0  =  10  (Test  6),  Uc/U0  = 
2.0  (Test  7),  and  a  reverse  direction  current,  U</U0  =  -2.0  (Test  8).  Again  the  model  was 


Table  4.2.  Model  hydrodynamic  input  for  wave  plus  current  sensitivity  study  of 
morphology  simulations.  _ _ _ _ 


Hydrodynamic  Input 

Test  #5 

Test#6 

Test  #7 

Test  #8 

U0  (m/s) 

0.35 

0.35 

0.35 

0.35 

T(s) 

4.79 

4.79 

4.79 

4.79 

Uc  (m/s) 

0.175 

0.35 

0.70 

-0.70 

A,  (m) 

0.23 

0.23 

0.23 

0.23 

V 

76 

76 

76 

76 

Field  Measure  Bedform 
(Year  day) 

315.9 

315.9 

315.9 

315.9 

run  for  twenty  wave  periods  for  each  test  case  with  the  results  shown  in  Figure  (4.3a-d). 
The  results  of  test  case  1  for  the  wave  only  forcing  are  included  in  the  figures  for 
comparison. 
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Figure  4.3.  DUNE2D  bedform  morphology  given  an  initial  bedform  (field  year  day 
315.9)  with  monochromatic  wave  (Uo=0.35m/s,  T=4.8s,  \j/=76 )  plus  current  forcing  after 
20  wave  periods.  (a)Test#5:Uc/Uo=0.5;  (b)  Test #6  Uc/Uo=1.0;  (c)  Test#7:  Uc/Uo=2.0;  (d) 
Test#8:  Uc/Uo=-2.0. 


For  Uc/U0  =  0.5  (Figure  4.3a),  the  current  has  little  effect  on  the  ripple  bed.  The 
20th  period  bedform  for  test  #5  is  very  similar  to  the  bedform  morphology  of  the  wave- 
only  test  case,  which  supports  Andersen’s  (1999)  model  results  for  wave  plus  low  current 
magnitudes  (IVU0  <  0.5)  having  minor  effect  on  bed  morphology.  For  Uc/U0  =  1.0 
(Figure  4.3b),  there  is  a  noticeable  migration  of  the  bed  in  the  direction  of  the  current. 
For  Uc/Uo  =  2.0  (Figure  4.3c),  more  sediment  is  removed  from  the  bed  and  placed  into 
suspension  and  greater  migration  of  the  bed.  Finally,  reversing  the  current  flow  direction 
180  degrees,  shifts  the  bed  migration  to  the  opposite  direction.  In  summary,  a  collinear 
current's  effect  on  the  bed  is  significant  at  values  where  the  current  magnitude  is  equal  to 
or  greater  than  the  wave  velocity  amplitude,  Uc/U0  >1.0. 

B.  MODEL  COMPARISON  WITH  SHOWEX  ’99  OBSERVATIONS 

ADUNE2D  morphology  simulation  of  the  SHOWEX  ’99  storm  event  described 
in  Chapter  III  was  conducted  to  assess  how  well  the  model  compares  with  real-world 
bedform  evolution  and  migration.  In  a  course  of  the  day,  field-measured  bedforms 
experienced  low  wave  energy  initially,  followed  by  large  wave  energy  generated  by  a 
coastal  storm  system  that  resulted  in  a  modified  bottom  bedform.  The  bedform  evolved 
from  a  mega-ripple  with  superimposed  small-scale  roughness  created  by  vortex  ripples  to 
a  flat,  smooth  bed  twenty-four  hours  later.  Refer  to  Figure  (3.7)  for  the  field  bedform 
morphology  time  series  produced  by  this  storm-forcing  event. . 

1.  Simulation  of  Storm  Event 

The  storm  event  is  simulated  by  two  separate  model  runs  representing  the  initial 

pre-storm,  low  wave  and  current  energy,  followed  by  the  conditions  at  the  peak  wave 

energy  plus  current  during  the  storm.  The  onshore,  cross-shore  current  was  low  for  both 
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simulations.  Most  of  the  current  was  oriented  in  the  along-shore  direction.  The 
orthogonal  currents  can  not  be  directly  included  into  the  simulation  since  DUNE2D  is 
limited  to  collinear  currents. 

An  approach  to  incorporate  the  effects  of  the  along-shore  currents  on  the  bedform 
is  to  use  an  effective  wave  velocity  amplitude  that  is  equivalent  to  the  combined  wave 
and  orthogonal  current  bed  shear  stress, 

'  ,  (4-1) 

=  Zt  + 

where  the  wave  bed  shear  stress  is  calculated  using  Equation  (2.14)  and  the  orthogonal 
current  stress  is  calculated  using  Soulsby’s  (1997)  equation  with  a  current  only  drag 
coefficient. 


^  =  fCDU2c 


(4.2) 


From  the  field  measurements,  the  calculated  along-shore,  current-only  shear  stress  is  1/3 
the  wave-only  bed  stress.  Therefore,  an  effective  stress  of  133%  of  the  wave-only  stress 
is  used  to  include  the  effects  of  the  orthogonal  current.  An  effective  wave  velocity 
amplitude,  U0,  is  obtained  from  this  effective  bed  stress  using  Equation  (2.14).  The 
results  of  the  modeled,  morphology  bedform  evolution  using  this  effective  wave  bed 
shear  stress  failed  because  the  model  generated  improper  equilibrium  ripple  lengths  and 
profiles.  Therefore,  the  along-shore  current  effects  on  the  bedform  are  assumed  not  to  be 
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important,  which  is  supported  by  laboratory  results  that  indicate  orthogonal  currents  do 
not  significantly  affect  maximum  shear  stresses  under  oscillatory  flow  (Soulsby  et  al, 
1993). 


Long-time  simulations  were  made  to  ensure  the  bed  would  reach  equilibrium, 
such  that  accurate  ripple  geometry  measurements  could  be  compared  to  the  field 
bedforms.  The  lower  energy  simulation  (Table  4.3.,  Test  #9,  v|/=16)  represents  the  wave 
conditions  at  the  beginning  of  the  storm,  near  year  day  315.9.  The  high-energy  waves 
conditions  on  year  day  316.4  (Table  4.3.,  Test#10,  \|/=341)  are  used  to  force  the  initial 
bedform  of  year  day  315.9  to  examine  the  morphology  evolution  and  migration  during 


the  peak  of  the  storm. 

Table  4.3.  Two  test  cases  in  storm  event  case  study  of  morphology  simulations. 


Hydrodynamic  Input 

Test  #9  (Low  Energy) 

Test#  10  (High  Energy) 

Uo  (m/s) 

0.16 

0.74 

T  (s) 

4.8 

9.6 

Uc  (m/s) 

0.11 

0.10 

V 

16 

341 

A,  (m) 

0.12 

1.13 

Bedform  (Year  day) 

315.9 

315.9 

Total  Simulation  Time  (s) 

9580 

9580 

The  simulated  bedform  morphology  time  series  are  shown  in  Figures  (4.4)  and 
(4.5a)  for  the  low  hydrodynamic  energy  with  the  field  bedform  morphology  evolution 
following  the  same  forcing  given  in  Figure  (4.5b)  for  comparison.  The  modeled 
morphology  for  the  high  energy  simulation  (test  #10)  is  illustrated  in  Figures  (4.6)  and 
(4.7a),  with  the  field  bedform  results  in  Figure  (4.7b)  for  comparison. 
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Figure  4.4.  Test#9:  DUNE2D  morphology  time  series  given  an  initial  bedform  (year 
day  315.9)  with  low  wave  energy  (Uo=0.16,  T=4.8,  Uc=-.ll,  \|/=16,  A=0.12m). 


Figure  4.5a.  Test  #9:  DUNE2D  morphology  bedform  snapshots  over  simulation  time 
resulting  in  a  vortex  ripple  field  at  equilibrium  state  (ripple  geometry:  A,Avg=0. 1 1  m,  p  Avg 
(height)  =  0.02m,  r)  Avg  A  Avg  (slope)  =  0.18,  XAvg/A=  0.92).  Migration:  4.2cm/hr  onshore. 


Figure  4.5b.  SHOWEX  ’99  Bedform  morphology  time  series  snap  shots  taken  from 
year  day  3 1 5.9  to  2.4  hours  later,  (superimposed  vortex  ripple  geometry:  XAvg=0.2m, 
PAvg  =  0.008m,  r\  Avg  /X  Avg  =  0.04)  Migration:  4.5cm/hr  onshore 
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day  315.9)  with  high  wave  energy  (Uo=0.74,  7=9.6,  Uc=-0.1,  Xj/— 341,  A=l.lm). 


Figure  4.7a.  Test  #10:  DUNE2D  Model  morphology  bedform  snapshots  over 
simulation  time  resulting  in  a  flat  mega-ripple  equilibrium  state  (ripple  geometry: 
^Average=l-8m,  T|  =  0.008m,  T|/A,  =  0.005,  7. Average/ A=  1.5). 


Figure  4.7b.  SHOWEX  ’99  Bedform  morphology  time  series  snap  shots  taken  from 
year  day  316.4(bedform  during  peak  wave  forcing)  to  3.6  hours  later.  (Ripple  geometry: 
7.Average=T-5m,  r|  =  0.02m,  r|/A,  =  0.02).  Migration;  2.2  cm/hr  on-shore. 
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2. 


Discussion  of  Storm  Event  Simulation 


The  low  energy  simulation  (test  #9)  reached  equilibrium  after  about  2000  wave 
periods  (+2.6hrs).  The  bedform  is  characterized  by  small-scaled  vortex  ripples  over  top  a 
larger  scale  ripple,  much  like  the  field  bedform  (Figure  4.5b).  The  average  wavelength  of 
the  vortex  ripples  at  equilibrium  state  is  0. 1 1  m,  with  an  average  height  of  0.02m.  A  ratio 
of  average  ripple  wavelength  to  orbital  displacement  is  approximately  equal  to  one.  This 
compares  qualitatively  well  with  the  observed  bedform  during  the  low  wave  energy  state 
and  with  Neilsen’s  (1981)  equation  (3.4.3)  for  V2 (k!  A),  given  v|/<200.  The  slope  of  the 
measured  vortex  ripples  superimposed  over  the  larger  scale  ripple  (Figure  4.5b)  is  0.04 
compared  with  the  modeled  ripples  slope  of  0.18.  Some  of  the  difference  is  due  to  the 
field  instrument’s  4cm  spatial  roll-off  in  the  horizontal  that  flattens  off  the  top  of  the 
small-scale  ripples  of  the  bed.  However,  since  the  observed  forcing  had  a  finite  spectral 
width  and  directional  spread,  the  SHOWEX  forcing  does  not  closely  match  the  single 
frequency  unidirectional  model  forcing.  The  migration  of  the  modeled  vortex  ripple  is 
4.2  cm/hr  onshore,  while  the  field  bedform  is  similar  at  4.5  cm/hr  onshore.  A  summary 
of  the  results  is  given  in  Table  (4.4). 


Table  4.4.  Storm  event  ripple  geometry  comparison  (DUNE2D  vs.  SHOWEX) 


Ripple  Geometry 

Output 

Low  Energy  Forcing 

High  Ener 

£?y  Forcing 

DUNE2D 

Test  #9 

SHOWEX 

(YR  315.9) 

DUNE2D 

Test  #10 

SHOWEX 

(YR  316.57) 

Xavz,  (m) 

0.11 

0.20 

1.8 

1.5 

riAve  (m) 

0.02 

0.008 

0.008 

0.02 

0.18 

0.04 

0.004 

0.01 

0.92 

1.7 

1.5 

1.4 

Migration  (cm/hr) 

4.2 

4.5 

2.2 
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The  high-energy  simulation  (test  #10)  reaches  equilibrium  in  about  1000  wave 
periods  (+2.6  hours),  faster  than  the  field  bedform  (+3.6  hours).  The  average  wavelength 
of  the  flat,  large-scale  ripple  is  1.8  m,  with  a  height  of  0.008  m.  The  ratio  of  ripple 
wavelength  to  orbital  displacement  is  around  1.5.  The  simulated  large-scale,  flat  ripple 
has  a  similar  dip  in  the  middle  of  the  bed  with  higher  peaks  towards  the  ends  of  the 
boundaries  as  the  field  bedform  (Figure  4.7b).  The  modeled  migration  rate  is  unresolved, 
while  the  measured  bed,  under  the  high  energy  forcing,  migrates  at  2.2  cm/hr  onshore. 


A  question  arises  with  the  ripple  geometry  for  the  high  energy  simulation  (Test 
#10)  as  to  whether  its  equilibrium  bedform  is  influenced  by  the  domain  size  of  the  grid, 
as  the  computed  equilibrium  ripple  wavelength,  ?iAvg,  appears  to  be  equal  to  the 
horizontal  domain  of  the  model  run.  (Figure  4.7a)  To  examine  the  question,  four  separate 
model  runs  were  performed  where  the  length  of  the  domain  of  the  initial  bedform  was  » 

varied,  as  summarized  in  Table  (4.5)  and  Figure  (4.8a-d). 


Table  4.5.  Grid  test  cases  to  verify  boundary  effects  on  morphology  simulations. 


Hydrodynamic  Input 

Test  #11 

Test#  12 

Test#13 

Test  #14 

Uo  (m/s),  T(s),  \j/ 

0.69  ,  9.6  ,  292 

A,  (m) 

1.05 

Initial  Bedform  (Year  day) 

315.9 

DUNE2D  Bedform  Horizontal 
Distance  Increase  from  Original  Bed 

30% 

50% 

100%  . 

300% 

Total  Simulation  Time  (s ) 

9580 

C  I  l  3  <  S 


Figure  4.8.  DUNE2D  initial  bedform  for  (a)Test  #1  l,(b)  Test  #12,  (c)  Test  #13  (d) 
Test  #14. 
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The  high  energy  model  simulations,  with  similar  forcing  as  test  #10,  produce 
different  equilibrium  ripple  wavelengths  than  measured  (Figure  4.9a-d).  The  model 


Figure  4.9.  DUNE2D  morphology  simulations  illustrating  bedform  evolution  at 
equilibrium  (+2.6hrs)  for  (a)Test  #ll(XAvg=1.8m),  (b)  Test  #12  (XAvg=2.2m),  (c)  Test  #13 
(A,Avg=l-5m).  (d)  Test#14  (unresolved  A,Avg ) 


output  for  the  30%  and  50%  bedform  enlargement  produce  single  mode  bedforms  with 
wavelengths  equal  to  the  total  horizontal  domain,  1.8m  and  2.2m,  respectively.  A 
doubled  domain  with  two  bedforms  produces  two  flatten  ripples  with  an  average 
wavelength  equal  to  1.5m,  about  half  of  the  total  horizontal  domain.  Test  #4  has  an 
unresolved  equilibrium  wavelength  due  to  the  large  perturbation  at  the  domain  ends.  The 
four  tests  show  the  equilibrium  wavelength  of  the  large-scale  ripple  that  evolves  under 
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high  energy  forcing  (forcing  capable  of  removing  vortex  ripples)  is  affected  by  the 
horizontal  domain  of  the  initial  bedform.  Therefore,  the  high  energy  results  can  be  only 
qualitatively  compared  with  the  field  bedform.  The  modeled  equilibrium  bedform  show 
similar  removal  of  vortex  ripples  and  a  flattened,  large-scale  ripple  for  high  energy 
forcing  as  the  field  measurements. 

C.  WAVE  SIMULATION  OF  FIELD  STORM  EVENT 

DUNE2D  is  limited  to  monochromatic  wave  forcing.  However,  the  field  data 
show  narrow-banded  wave  forcing  with  some  groupiness  present.  To  simulate  the 
groupiness  of  the  field  measured  waves,  a  sequence  of  monochromatic  wave  model  runs 
was  performed.  Each  successive  model  run  was  represented  by  the  equivalent  wave 
velocity  amplitude  and  period  of  the  measured  waves.  The  duration  of  each  run  was 
equal  to  the  measured  duration  of  the  section  of  real  wave  condition,  so  that  a  piecewise 
narrow  band  forcing  time  series  could  be  constructed.  Additionally,  each  model  run  was 
initiated  by  the  previous  final  morphology  bedform  output  and  flow  solutions. 

However,  the  method  failed  when  trying  to  initialize  the  successive  runs  with  the 
previous  bedform  output  and  flow  solutions.  For  the  case  of  a  moveable  bed  as  used  in 
these  simulations,  the  initial  bedform  is  normalized  by  a  scaling  parameter,  D,  which  is 
based  on  the  monochromatic  forcing  input  (U0  and  T).  The  model  generates  a  grid  from 
this  initial,  normalized  bedform.  When  the  successive  model  runs  are  initiated  with  the 
previous  output  bedform  and  flow  solutions,  the  model  is  unable  to  resolve  the  new  flow 
solutions  due  to  different  scaling  values  used  to  generate  dissimilar  grids  between  the  two 
model  runs.  Further  work  is  required  to  incorporate  the  amplitude  changes  associated 

with  groupiness  waves  that  are  representative  of  the  real  world. 

40 


V.  CONCLUSION 


The  migration  and  evolution  of  small-scale  bedforms  was  simulated  using 
DUNE2D  based  on  field  hydrodynamic  forcing  of  the  bottom  boundary  layer.  The  model 
simulates  flow,  sediment  transport,  and  morphology  in  the  bottom  boundary  layer  to 
enable  bedform  comparisons  with  the  field  data.  The  initial  cross-shore  bedform  used  in 
the  model  closely  matched  the  measured  SHOWEX  acoustic  altimeter  data.  However,  it 
is  required  to  manipulate  the  measured  bedform  to  conform  to  the  cyclic  boundary 
requirement.  At  low  energy  forcing,  the  model  generated  vortex  ripples.  The 
equilibrium  wavelengths  of  the  rippled  bedforms  were  near  the  orbital  diameter  for  the 
oscillatory  wave  forcing.  SHOWEX  bedform  changes  under  low  wave  plus  collinear 
current  conditions  resulted  in  minor  changes  of  the  vortex  ripple  fields.  Bedform 
migration  rates  of  the  model  were  similar  to  the  field  migration  rates. 

DUNE2D  simulations  of  high  energy  forcing  over  a  ripple  bed  resulted  in  greater 
change  in  the  bed  over  a  shorter  duration  of  time  compared  with  the  lower  energy  forcing 
simulation.  Like  the  field  data,  the  modeled  data  under  strong  forcing  removed  smaller 
scale  vortex  ripples  and  redistributed  the  sediment  into  a  larger  scale  ripple  with  a  large 
portion  of  sediments  in  suspension  above  the  bed.  Final  equilibrium  bedforms  consisted 
of  larger  scale,  flattened  ripples  with  wavelengths  reaching  beyond  orbital  excursion 
values.  However,  the  horizontal  domain  of  the  initial  bedform  and  the  number  of  initial 
low-mode  ripples  was  found  to  be  somewhat  influenced  by  the  equilibrium  ripple 
wavelength  of  the  large-scale  ripple  that  evolved  under  high  energy  forcing. 
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In  nature,  sediment  sizes  and  distribution  varies  spatially  and  temporally.  The 
model  simulations  are  limited  to  a  single  sediment  grain  diameter  designation,  which  in 
turn  affects  the  Shield  parameter  calculations.  Ultimately,  the  Shield’s  parameter 
calculations  relate  to  sediment  transport  and  bedform  morphology  accuracy. 

The  oceanic  wave  forcing  is  not  monochromatic  and  unidirectional,  so  detailed 
comparisons  with  DUNE2D  model  are  not  possible.  However,  the  modeled  bedform  data 
was  compared  qualitatively  with  changes  in  the  measured  bedforms  (Table  4.4). 
Differences  can  be  ascribed  to  model  input  limitations,  which  include  simple  sinusoidal 
wave  forcing,  collinear  current  directionality,  periodic  bedform  boundaries,  and  single 
sediment  grain  diameter  selection.  Efforts  focused  on  generalizing  the  wave  forcing  to 
.simulate  actual  field  conditions  and  initial  bed  configuration.  Simulation  of  the  storm 
event  through  sinusoidal  flow  over  bedforms  restricted  actual  simulation  of  real  waves. 
The  higher  energy  SHOWEX  velocity  fields  were  characterized  by  groupy  waves.  Wave 
simulation  could  not  duplicate  the  groupiness  observed  in  the  field.  The  inability  to  vary 
current  direction  within  the  model  reduced  proper  sediment  transport.  During  the  storm 
event,  strong  along-shore  current  existed.  Simulation  of  collinear  currents  demonstrated 
substantial  transport  when  currents  are  greater  than  the  velocity  amplitude,  Uc/U0  >  1.0. 
Despite  the  inability  to  establish  arbitrary  wave  forcing  and  directional  currents,  the 
model  bedform  tendencies  based  on  the  simple  forcing  provided  good  qualitative 
agreement  with  field  bedforms.  Continued  work  on  the  capability  of  the  model  to 
incorporate  real  world  wave  and  current  forcing  is  encouraged.  The  model  provides  a 
good  basis  for  which  to  better  understand  the  probable  bed  migration  and  evolution  over 
time  given  near  monochromatic  wave  forcing  in  the  nearshore. 


42 


APPENDIX  A.  DUNE2D  MODEL  SOLUTION 


This  appendix  provides  detailed  discussion  of  DUNE2D  model  solutions.  The 
following  sections  are  divided  between  model  input  and  output  descriptions. 

A.  MODEL  INPUT 


Hydrodynamic  forcing,  geophysical  properties  of  the  sediment,  numerical 
schemes,  grid  definitions,  boundary  conditions,  transport  methods,  and  morphology  setup 
are  discussed  below. 

1.  Hydrodynamic  and  Geophysical  Input 


The  model  is  forced  by  monochromatic  sinusoidal  wave  motion.  Wave  motion  is 

described  by  the  amplitude  of  the  oscillatory  wave  velocity,  U0,  and  period,  T. 

Implementation  of  a  current  in  the  same  direction  as  the  wave  is  an  option.  The  scaling 

depth,  D,  is  applied  to  most  input  terms  to  provide  non-dimensionality  where 

D  =  2m  =  U0T  (Al) 

The  model  is  based  on  turbulent  boundary  flow.  The  flow  is  characterized  by  the 

Reynolds  number  computation  as  described  below: 


1C.M, 


v  =  kinematic  viscosity 


(A.2) 


V 

Sediment  description  is  limited  to  a  homogeneous  bottom  type,  described  by  the 
medium  grain  diameter,  dso.  The  grain  diameter  is  incorporated  into  the  model  through 
the  roughness  of  the  bottom,  kn  (Equation  2.10).  (Fredsoe  and  Deigaard,  1992)  Porosity 
of  the  sediment  is  set  at  0.4.  The  ratio  of  sediment  density,  ps,  to  water  density,  p,  is 
constant  at  2.65,  based  on  in-situ  sand  samples. 
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Model  simulation  time  is  based  on  total  number  of  wave  periods.  Time  steps,  At, 
are  a  function  of  the  number  of  iterations  for  each  wave  period.  Both  time  functions  are 
described  by  the  DUNE2D  model  input  variable  WAVE. 


Table  A.  1 .  Hydrodynamic  and  geophysical  input. 


DUNE2D 

Variable 

Physical  Parameter 

ROUG 

Non-dimensional  Nikuradse  constant,  kn/D 

ANGU 

Non-dimensional  depth,  D/a 

WAVE 

Total  number  of  wave  periods  &  iterations  per  period 

TURB 

Reynolds  number,  Re 

2.  Numerics  and  Scheme  Setup 

A  finite  volume  numerical  method  is  used  to  discritize  the  general  flow  equations 
in  section  (2.3).  The  flow  is  assumed  to  be  incompressible.  The  pressure  term  in  the 
Reynolds-averaged  Navier-Stokes  (RNS)  equations  is  calculated  using  the  PISO 
algorithm  (Patankar,  1980).  The  PISO  algorithm  is  an  iterative  process  used  to  determine 
the  pressure  values  at  each  grid  point.  The  user  provides  the  number  of  times  the  PISO 
algorithm  is  implemented.  For  flows  that  are  difficult  to  converge,  a  higher  number  of 
iterations  are  required.  Andersen  emphasizes  the  difficulty  of  solving  the  pressure  term 
within  the  RNS  equations. 

Zijlema’s  (1996)  ISNAS  scheme,  a  high  order  spatial  discretization  method,  is 
used  to  remove  numerical  diffusion  errors.  Implicit  time  discretization  methods  are 
utilized,  except  for  the  advection  terms  where  semi-implicit  techniques  are  implemented. 


44 


Table  A. 2.  Numerical  Input. 


DUNE2D 

Variable 

Physical  Parameter 

SCHM 

Discretization  using  ISNAS  scheme 

PISO 

PISO  algorithm  parameter 

TURB 

k-g)  turbulent  closure  model 

3.  Grid  Input 

For  morphology  simulations,  the  grid  is  generated  using  a  transfinite  interpolation 
technique.  The  transfinite  interpolation  technique  produces  finer  resolution  at  the  bed, 
while  vertically  stretching  the  grid  points  out  as  the  distance  from  the  bed  increases.  A 
specific  constant,  GMET,  is  used  to  define  the  degree  of  stretching.  Large  GMET  values 
correspond  to  increased  stretching.  Additionally,  the  total  number  of  vertical  coordinate 
points,  M,  is  a  function  of  the  transfinite  interpolation  method. 

The  user  can  specify  bedform  through  an  ASCII  file  or  have  DUNE2D  generate  a 
bottom  using  several  predefined  bedform  shapes.  For  user-specified  bedforms,  the 
ASCII  file  contains  information  about  the  number  of  horizontal,  N,  and  vertical,  M,  grid 
points.  Also,  the  file  contains  the  bedform  grid  values  normalized  by  the  scaling  term,  D. 
An  example  of  a  grid  file,  test.grd,  is  given  in  Appendix  (B). 

The  resolution  of  the  grid  spacing  is  set  by  the  bed  roughness,  which  is  a  function 
of  the  bedform  profile,  Nikuradse  constant,  sediment  grain  diameter  and  wave  forcing 
input.  From  Fredsoe  and  Deigaard,  the  relationship  used  to  determine  the  resolution 
parameter,  RESO,  requires  the  calculation  of  the  maximum  frictional  velocity,  u  , 
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u*  -  0.1-j2U„\  — 


.  2.5  u 

RESO  =  integer  — d50  — 
0.4  50  v 


(A- 3) 


RESO  is  used  in  the  grid  generating  procedure.  Large  RESO  values  are  indicative 


of  the  higher  grid  resolution. 


DUNE2D 

Variable 


GRID 


RESO 


GMET 


Table  A.3.  Grid  Input. 


Physical  Parameter 


N:  horizontal  total  grid  pts;  M:  vertical  grid  pts 


Grid  Resolution  parameter 


TFI  grid  generator  stretching  parameter 


An  example  of  a  generated  grid  using  the  transfmite  interpolation  technique  over 
a  user-defined  bedform  is  shown  in  Figure  (A.  1). 


Morphology  Grid  File  -  Simple  Transfinite  Interpolation  Grid  Generator 


ifiipf 


m  i 


GMETr=L2980 


0  0.05  0.1  0.15  0.2  0.25  0.3  o:35 

x/D  (m) 

Figure  A.  1 .  DUNE2D  grid  using  the  transfinite  interpolation  method. 
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4. 


Grid  Sensitivity  Study 


A  grid  sensitivity  test  of  the  flow  performed  by  Andersen  (1999)  is  summarized  in 
Table  (A.6).  Four  parameters  were  tested:  y+,  M,  Ax^st/a,  and  D.  The  first  term,  y+,  is 
the  distance  from  the  bottom  boundary  to  the  first  grid  point.  Flow  around  the  crest  of 
the  ripples  was  found  to  be  sensitive  to  the  bed  slope.  A  stretched  grid  was  applied  to 
enhance  the  number  of  points  in  that  region.  Therefore,  the  value  Axcrest,  the  horizontal 
spacing  between  grid  points  around  the  crest,  was  tested.  D  is  the  normalized  depth  of 
the  flow.  The  tests  provide  guidance  on  model  performance,  in  particular  on  flow 
convergence  and  numerical  stability. 


Table  A.4.  Grid  Sensitivity  of  the  Flow  Test  Results. 


+ 

y 

<=0.40 

NxM 

40X30 

AXcrest/ & 

<=0.012 

D 

>=5.0*a 

Andersen  also  performed  a  grid  sensitivity  test  for  the  suspended  load.  By 
varying  the  total  vertical  grid  points  between  40,  60,  and  100,  little  was  achieved  in  grid 
convergence  of  the  suspended  concentration  above  the  bed. 

5.  Model  Boundary  Conditions 


The  lateral  boundary  conditions  are  periodic.  Bottom  boundary  conditions  are: 


o(  ) 

— —  =  0  u,  w  =  0  at  the  z  =  D 
8z 
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The  no-slip  condition  creates  a  viscous  sub-layer.  The  top  boundary  is  described  by  a 

symmetric  boundary  condition. 

6.  Sediment  Transport  Related  Setup 

Bedload  and  bed  concentration  formula  are  based  on  Engelund  and  Fredsoe 
formulation,  which  is  described  in  Chapter  II.  Angle  of  repose  of  the  sediment  is 
constant  at  0.65.  Critical  Shields  parameter  is  0.05  for  all  model  runs. 

The  sediment  fall  velocity  is  calculated  using  Fredsoe  and  Deigaard  s  empirical 
relationship  for  sand  using  an  iterative  process  of  Equation  (2.17). 

Suspended  load  along  the  lateral  boundaries  is  treated  as  a  true  cyclic  boundary. 
Additionally,  wall  boundary  techniques  are  used  to  separate  the  suspended  load  from 
bedload.  The  bed  concentration  is  assumed  to  be  based  on  normal  diffusion. 

Table  A.  5 .  Sediment  Transport  Input. 

Physical  Parameter 

Bedload  and  Bed  concentration  formula  using  Engelund-Fredsoe 
_ formula  _ _ _ 

_ Non-dimensional  settling  velocity  of  the  sediment _ 

_ Wall  boundary  option  for  suspended  load _ 

_ Normal  diffusion  of  bed  concentration _ _ 

_ True-cyclic  suspended  load _ _ _ 

7.  Morphology  Setup 

Once  the  morphology  mode  is  turned  on,  only  a  few  additional  DUNE2D 
variables  are  required.  The  QUICK  scheme  with  top  smoothing  is  applied  to  the  solution 


DUNE2D 
Variable 
BEDL,  BEDC 

FALL 

SSLO 

SMIX 

SSCY 


of  the  continuity  equation  in  the  morphology  section  below.  This  is  third  order  accurate. 
(Andersen,  1999).  The  smoothing  is  based  on  a  running  average. 


Table  A.6.  Morphology  Input. 


DUNE2D 

Variable 

laUiU  TV.V.  ; - | - - — - 

Physical  Parameter 

MSCH 

Leonard  (T979Vs  Quick  Scheme  for  continuity  equation 

OSMO 

Smoothing  is  iterated  4  times  over  bed  and  at  ripple  crests 

B.  MODEL  OUTPUT 

Figure  (A.2)  summarizes  output  parameters  from  each  of  the  three  modules. 


Figure  A.2.  DUNE2D  Module  Output  Diagram. 


The  morphology  mode  simulation  produces  time-averaged  quantities  over  a  single 
wave  period.  Table  (A. 7)  summarizes  the  different  output  variables  that  are  saved  to 
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ASCn  files  by  DUNE2D.  Most  of  the  output  is  non-dimensional.  Transformation  of 
these  quantities  into  meaningful  dimensional  output  will  be  discussed  in  Appendix  (C). 


Table  A.7.  Morphology  Mode  Output.  (Period  Average  Mean). 


A  auJLW  /  . 

Output 

Physical  Parameter 

<Tb>t 

Bed  shear, time  averaged  over  one  wave  period 

<0>t 

Shield  parameter,  time  averaged  over  one  wave  period 

<qb>t 

Bedload  flux,  time  averaged  over  one  wave  period 

<Qs>t 

Suspended  load  flux,  time  averaged  over  one  wave  period 

<qt>t 

Total  load  flux,  time  averaged  over  one  wave  period 

_ _  XL  .  - - - - 

<4>t>t 

Non-dimensional  total  load  flux,  time  averaged  over  one  wave 

period 

<C>t 

Bed  concentration,  time  averaged  over  one  wave  period 

<u>t 

Horizontal  velocity,  time  averaged  over  one  wave  period 

<W>t 

Vertical  velocity,  time  averaged  over  one  wave  period 

<Dh/dt>t 

Change  in  bed  form  w/  time,  time  averaged  over  one  wave  period 

<h>t 

Bed  form  height,  time  averaged  over  one  wave  period 
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APPENDIX  B.  DUNE2D  MODEL  APPLICATION 


This  appendix  provides  further  information  on  the  application  of  DUNE2D  code. 
The  first  section  describes  general  DUNE2D  Fortran  file  information,  as  well  as  model 
initiation.  The  following  section  emphasizes  input  file  examples. 

A.  FORTRAN  CODE  FILES 


DUTME2D  is  a  modular  Fortran  77  code  program  consisting  of  (16)  *  f  programs, 
(9)  *  i  programs,  and  a  makefile.  The  code  is  divided  into  logical  modules  designed 
around  the  model  input,  flow  calculations,  bedload  transport  computations,  morphology 
solutions,  and  output.  The  makefile  is  designed  to  compile  all  the  programs  into  a  single 
executable  program  called  dune.  A  list  of  programs  is  given  by  Table  (B.  1). 


In  the  Unix  environment,  typing  dune  followed  by  a  string  four-letter  project 
name,  e.  g.  “test”,  will  execute  the  simulation  application.  Other  files  are  required  to  run 
DUNE2D.  These  will  be  discussed  in  the  next  section. 

Table  B.  1 .  Fortran  DUNE2D  Files.  {Program  areas:  a)I/0:  input/output  b)F:  flow 

.  4  •  _ _ l  JMV  /f  •  /-'v  t-^'a  V\  r\rr\r  MPt/R-  oriH/hnimdarvl 


seaimenx  iranspun  u 

b  .f :  G/B 

cmd.f :  I/O 

O -  -  -J 

dune2d.f :  main 

f.f :  F 

g.f.G/B 

if:  I/O 

io.f :  I/O 

isolines.f :  G/B 

k.f :  F 

ko.f :  F 

m.f :  F 

profile.f :  F 

p.f :  F 

r.f :  F 

ripple,  f :  G/B 

susp.f :  S 

s.f :  S 

t.f :  F 

wave.f :  F 

wavemorph.f :  M 

51 


B.  MODEL  INITIATION  FILES 

Once  an  error-free  compiled  version  of  dune  exists,  simulation  can  begin  with  the 
proper  input  files  in  place.  The  files  have  *.inp,  *.grd,  and  *  fort  file  extension,  which 
represent  input,  grid,  and  format  ASCII  files.  The  three  files  must  have  the  same  four 
letter  project  name  that  is  typed  after  dune,  e.  g.  test.inp,  test.grd,  and  test.fmt. 
Additionally,  special  input  files,  Special.inp  and  Equations. inp,  must  be  included  in  the 
simulation  folder.  These  two  files  are  used  to  initite  DUNE2D  simulation 

1.  Input  File 

The  input  file,  test.inp,  describes  the  necessary  input  variables  used  to  describe 
model  calculation,  initiation,  hydrodynamic  forcing,  total  time  of  integration,  and  other 
requirements  used  in  flow,  bedload  transport,  and  morphology  modules  described  in 
Chapter  II.  Appendix  (A)  provides  greater  detail  on  the  specific  requirements  of  the 
input  file.  An  example  is  provided  at  the  end  of  this  section. 

2.  Grid  File 

The  grid  file,  test.grd  or  test.xyb,  is  a  two-column,  ASCII  file  that  describes  the 
initial  bedform  grid  point  locations.  User-defined  grids  are  described  by  *.xyb  files  that 
are  the  same  format  as  *.grd  files.  An  example  has  been  included. 

3.  Format  File 

Finally,  the  format  file  is  used  to  describe  the  boundary  conditions  of  the  model. 
Since  all  model  runs  for  this  thesis  and  the  work  by  Andersen  (1999)  only  use  a  single 
block,  the  format  file  was  not  required  in  the  model  simulation.  For  more  information  on 
the  model  boundary  conditions  refer  to  Appendix  (A). 
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Fxamole  of  input,  file:  testing 

#  General  commands:  test.inp  03-Nov-2000 
READ,  0,  0.0 

COUT,  0,  0.0 
GOUT,  0,  0.0 
SCRE,  10,  0.0 
PISO,  3,  0.0 

#  Commands  for  the  solver: 

ZTOP ,  1,  0.0 

TURB,  4,  5.87e+05 
XNIT,  2,  0.0 
SCHM,  4,  1.0000 
SKEW,  1,  0.0 

#  Setup  of  the  grid: 

COOR,  1,  0.0 

GRID,  100,  40.0000 
GMET,  4,  1.2980 
RESO,  13,  0.0 

#  Bed  load  transport: 

BEDL,  1,  0.0 

ROUG,  0,  0.0001 
TTOP ,  1,  0.0863 
SPHI,  0,  0.5770 

#  Suspended  transport: 

SUSP,  1,  0.0 

SSCY ,  3,  0.0 
SSLO,  1,  0.0 
BEDC,  0,  0.0 
SMIX,  0,  0.0 
FALL,  2,  0.0168 

#  Setup  of  the  wave: 

DYNA,  1,  0.0 

ANGU,  0,  6.2832 
WAVE,  20,  1000.0000 
WEND,  0,  1.00e-24 
MEAN,  1,  0.0 

#  Morphology  Module: 

MORP,  1,  0.0 

MOUT,  150,  0.0 
MOIT ,  5 ,  0.0 
MSCH,  8,  0.0 
QSMO,  4,  4.0000 
MRLX,  0,  1.0000 

#  User  defined  variables: 

PROB,  0,  0.0 
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120  40 

0.00000000  -0.00265241 

0.00261905  -0.00285512 

0.00523809  -0.00189552 

0.00785714  -0.00112030 

0.01047619  -0.00080800 

0.01309524  -0.00013576 

0.01571429  0.00004675 

0.01833333  0.00102395 

0.02095238  0.00214214 

0.02357143  0.00334799 

0.02619048  0.00452660 

0.02880952  0.00567079 

0.03142857  0.00685681 

0.03404762  0.00804007 

0.03666667  0.00920828 

0.03928571  0.01034211 

0.04190476  0.01155320 

0.04452381  0.01285784 

0.04714286  0.01420774 

0.04976190  0.01547332 

0.05238095  0.01666995 

0.05500000  0.01779334 

0.05761905  0.01893293 

0.06023810  0.02016495 

0.06285714  0.02150976 

0.06547619  0.02288504 

0.06809524  0.02416856 

0.07071429  0.02539719 

0.07333333  0.02664316 

0.07595238  0.02788491 

0.07857143  0.02913498 

0.08119048  0.03033831 

0.08380952  0.03156390 

0.08642857  0.03283032 

0.08904762  0.03408979 

0.09166667  0.03537232 

0.09428571  0.03657105 

0.09690476  0.03779511 

0.09952381  0.03902313 

0.10214290  0.04032772 

0.10476190  0.04169318 

0.10738100  0.  _ 


Program  break.  Note  the  first  line  is  the  total  number  of  horizontal  and  vertical 


points. 

4.  Summary 

Summaries  of  the  required  files  for  simulation  are  given  by  Table  (B.2).  These 

files  should  all  be  placed  in  the  same  Unix  folder. 
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APPENDIX  C.  DUNE2D  INPUT/OUTPUT  TRANSFORMATION 

This  appendix  provides  specific  input  and  output  procedures  used  to  implement 
DUN0E2D  code.  An  overview  of  the  process  by  which  data  is  entered  into  the  DUNE2D 
Fortran  format  input  files  as  well  as  transformation  of  the  output  files  into  meaningful, 
dimensional,  as  well  as,  graphical  results  is  shown  in  Figure  (C.  1). 


Figure  C.l.  General  input/output  flow  diagram  of  modeling  simulation. 

A.  GENERAL  I/O  PROCEDURES 

The  first  step  to  a  morphology  simulation  focuses  on  determining  the  wave  and 
current  forcing  dynamics.  A  Matlab  program  gathers  user  defined  input  concerning  the 
ocean  dynamics  and  sediment  mean  grain  diameter.  An  *.inp  file  is  created.  Second,  a 
bedform  *.grd  file  is  generated  from  field  or  experimental  bedform  data.  A  schematic 
diagram  of  the  MATLAB  programs  used  to  generate  the  input  file  data  is  shown  in 

Figure  (C.2). 

57 


Figure  C.2.  DUNE2D  input  data  flow  of  MATLAB  programs. 

Once  all  the  input  files  are  created.  The  first  step  of  the  morphology  simulation  is 
to  run  a  hot  file  simulation  run.  The  main  purpose  of  the  hot  file  run  is  to  generate  an 
initial  solution  to  the  flow  field  over  a  flat  bed.  The  output  generated  by  DUNE2D 
includes  a  *.hot  file  at  the  end  of  each  wave  period  of  simulation.  Normally,  the  hot  file 
is  generated  over  four  complete  wave  periods.  The  *  hot  output  file  is  used  by  the  actual 
morphology  simulation  to  help  initiate  the  model  flow  solutions.  If  the  hot  file  is  not 
included,  the  morphology  run  from  a  “cold”  start  will  typically  not  have  convergence  of 

the  flow  solution. 

During  the  model  simulation,  output  in  ASCII  format  of  specific  single  wave 

period  averaged  parameters  (Table  A.7)  are  placed  into  *.meanbed  files  in  the  Unix 

simulation  folder.  MATLAB  programs  were  developed  to  produce  graphical  data  output. 

A  summary  of  the  programs  and  output  options  are  included  in  Figure  (C.2). 
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Dune2d  grid  generator  for  initial 
bedform 

Bedform  morphology  snapshots  at  the  end 
of  each  wave  period  output  (bedform  timeseries) 

mb  Bedform  morphology  movie  made  from  mean  bed 
ESiH  snapshots  at  each  wave  period  (bedform  timeseries) 


S|gflll|  Input  data  (Um,  T,  d50)  comparisons  based  on 
Mlill  bedform  morphology  snapshots  at  1st  and  20th  v 
period  for  each  different  model  run. 

Initial  bedform  used  as  input  to  Dune2d 


Bedform  morphology  over  X  number  of  wave  periods,  as 
well  as.  bedload  (qb).  suspend  load  (qs).  total  load  (qt) 
fluxes,  bedshear  (tauj  for  the  1st, 5*  10th,  and  20*  wave  period 

Bedform  morphology  over  X  number  of 
wave  periods  (bedform  waterfall  plots) 


Ficrure  C.2.  DUNE2D  output  data  flow  of  MATLAB  programs. 


59 


B.  FIELD  DATA  TRANSFORMATION 

More  detail  on  the  specific  MATLAB  programs  used  to  transform  field  data  into 
DUNE2D  input  variables  and  files  are  provided  in  Figure  (C.3). 


Figure  C.3. 


MATLAB  programs  used  to  transform  field  data  into  model  input. 


X.  Bedform  Transformation 

MATLAB  program  make  xyb.m  produces  non-dimensional  bedform  (Figure 
C.4a-f.)  in  the  form  of  an  ASCII  file  *.xyb,  as  seen  in  Appendix  (B).  Initially  the  field 
bedform  data  points  (Figure  C.4a.)  are  re-oriented  so  that  the  left  end  point  starts  at  0 
meters,  in  both  height  and  horizontal  distance  (Figure  C.4b.).  The  boundary  conditions 
require  the  bedform  to  be  periodic  on  each  end.  This  requires  ficticous  points  to  be  added 
to  the  right  end  of  the  shifted  bedform.  The  bedform  was  extended  an  additional  30%  of 
the  original  bed  to  produce  a  gentler  slope  up  towards  the  zero  height  line.  Twenty 
additional  points  were  used  to  linearly  interpolate  from  the  last  original  data  point  on  the 
right  side  up  to  the  height  of  the  point  on  the  left  boundary  (Figure  C.4c).  To  ensure 
smooth  ends  for  the  fully  periodic  boundaries,  a  10%  cosine  filter  is  applied  to  the 
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extended  bedform  (Figure  C.4d).  The  filter,  when  applied,  maintains  most  of  the  original 
roughness  of  the  bedform  except  at  the  ends,  where  the  smoothing  takes  place  (Figure 
C.4e).  The  final  step  is  to  make  the  length  and  height  of  the  bed  non-dimensional  by 
dividing  through  by  a  length  scale,  D  (Figure  C.4f).  Further  information  on  calculation 

of  D  can  be  found  in  Appendix  (A). 


Figure  C.4.  Transformation  of  cross-shore  field  bedform  data  to  DUNE2D  non- 
dimensional  bedform  data. 
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2. 


Monochromatic  Wave  Dynamic  Forcing  Transformation 


MATLAB  program  waveforce.m  is  used  to  determine  the  monochromatic 
sinusoidal  wave  forcing  described  by  the  oscillatory  velocity  amplitude  and  peak  wave 
period  that  are  incorporated  into  the  input  file,  *.r«p.  From  the  field  measured  it  (along¬ 
shore)  and  v  (cross-shore)  velocity  profiles  of  the  BCDV  profiler,  the  power  spectrum  can 
be  found  to  determine  the  peak  period,  T,  of  the  wave  forctng.  The  oscillatory  wave 
velocity  amplitude  is  calculated  as  the  equivalent  sinusoidal  velocity  amplitude  that  gives 

the  same  wave  variance  as  measured, 

_  (C.l) 

U0  =  J2^ct;+ct; 

where  ou2  and  cv2  are  the  along-shore  and  cross-shore  velocity  component  variances. 
This  calculation  is  performed  to  get  the  best  directional  representation  of  the  wave 
forcing  over  the  bed. 

Finally,  mean  values  of  the  u  and  v  velocity  time  series  represent  the  along-shore 
and  cross-shore  currents,  respectively. 
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